%pylab inline
colors = [[0.1,0.6,0.5], [0.9,0.2,0.5], [0.9,0.5,0.1]]
Populating the interactive namespace from numpy and matplotlib
from PyDynamic.uncertainty.propagate_DFT import GUM_DFT,GUM_iDFT
from PyDynamic.uncertainty.propagate_DFT import DFT_deconv, AmpPhase2DFT
from PyDynamic.uncertainty.propagate_DFT import DFT_multiply
#%% reference data
ref_file = np.loadtxt("DFTdeconv reference_signal.dat")
time = ref_file[:,0]
ref_data = ref_file[:,1]
Ts = 2e-9
N = len(time)
#%% hydrophone calibration data
calib = np.loadtxt("DFTdeconv calibration.dat")
f = calib[:,0]
FR = calib[:,1]*np.exp(1j*calib[:,3])
Nf = 2*(len(f)-1)
uAmp = calib[:,2]
uPhas= calib[:,4]
UAP = np.r_[uAmp,uPhas*np.pi/180]**2
#%% measured hydrophone output signal
meas = np.loadtxt("DFTdeconv measured_signal.dat")
y = meas[:,1]
# assumed noise std
noise_std = 4e-4
Uy = noise_std**2
Consider knowledge about the measurement system is available in terms of its frequency response with uncertainties associated with amplitude and phase values.
\mathbf{H} = (\vert H(f_1) \vert, \ldots, \angle H(f_N))
u_H = (u_{\vert H(f_1) \vert}, \ldots, u_{\angle H(f_N)})
figure(figsize=(16,8))
errorbar(f * 1e-6, abs(FR), 2 * sqrt(UAP[:len(UAP) // 2]), fmt= ".-", alpha=0.2, color=colors[0])
xlim(0.5, 80)
ylim(0.04, 0.24)
xlabel("frequency / MHz", fontsize=22); tick_params(which= "both", labelsize=18)
ylabel("amplitude / V/MPa", fontsize=22);
figure(figsize=(16,8))
errorbar(f * 1e-6, unwrap(angle(FR)) * pi / 180, 2 * UAP[len(UAP) // 2:], fmt= ".-", alpha=0.2, color=colors[0])
xlim(0.5, 80)
ylim(-0.2, 0.3)
xlabel("frequency / MHz", fontsize=22); tick_params(which= "both", labelsize=18)
ylim(-1,1)
ylabel("phase / rad", fontsize=22);
The measurand is the input signal \mathbf{x}=(x_1,\ldots,x_M) to the measurement system with corresponding measurement model given by
y[n] = (h\ast x)[n] + \varepsilon[n]
Input estimation is here to be considered in the Fourier domain.
The estimation model equation is thus given by
\hat{x} = \mathcal{F}^{-1}\left( \frac{Y(f)}{H(f)}H_L(f) \right)
with - Y(f) the DFT of the measured system output signal - H_L(f) the frequency response of a low-pass filter
Estimation steps
- DFT of y and propagation of uncertainties to the frequency domain
- Propagation of uncertainties associated with amplitude and phase of system to corr. real and imaginary parts
- Division in the frequency domain and propagation of uncertainties
- Multiplication with low-pass filter and propagation of uncertainties
- Inverse DFT and propagation of uncertainties to the time domain
With the DFT defined as
Y_k = \sum_{n=0}^{N-1} y_n \exp(-\mathrm{j} k\beta_n)
with \beta_n = 2\pi n /N, the uncertainty associated with the DFT outcome represented in terms of real and imaginary parts, is given by
U_{Y} = \left( \begin{array}{cc} C_{\cos} U_y C_{\cos}^{\sf T} & C_{\cos} U_y C_{\sin}^{\sf T} \\ (C_{\cos} U_y C_{\sin}^{\sf T})^{\sf T} & C_{\sin} U_y C_{\sin}^{\sf T} \end{array}\right)
Y,UY = GUM_DFT(y,Uy,N=Nf)
figure(figsize=(18,6))
subplot(121)
errorbar(time*1e6, y, sqrt(Uy)*ones_like(y),fmt=".-")
xlabel("time / µs",fontsize=20); ylabel("pressure / Bar",fontsize=20)
subplot(122)
errorbar(f*1e-6, Y[:len(f)],sqrt(UY[:len(f)]),label="real part")
errorbar(f*1e-6, Y[len(f):],sqrt(UY[len(f):]),label="imaginary part")
legend()
xlabel("frequency / MHz",fontsize=20); ylabel("amplitude / au",fontsize=20);
In practice, the frequency response of the measurement system is characterised in terms of its amplitude and phase values at a certain set of frequencies. GUM uncertainty evaluation, however, requires a representation by real and imaginary parts.
H_k = A_k \cos(P_k) + \mathrm{j} A_k\sin(P_k)
GUM uncertainty propagation
C_{RI} = \left( \begin{array}{cc} R_{A} & R_{P} \\ I_{A} & I_{P} \end{array}\right) .
U_H = C_{RI} \left( \begin{array}{cc} U_{AA} & U_{AP} \\ U_{AP}^{\sf T} & U_{PP} \end{array} \right) C_{RI}^{\sf T} = \left( \begin{array}{cc} U_{11} & U_{12} \\ U_{21}^{\sf T} & U_{22} \end{array} \right) .
H, UH = AmpPhase2DFT(np.abs(FR),np.angle(FR),UAP)
Nf = len(f)
figure(figsize=(18,6))
subplot(121)
errorbar(f*1e-6, H[:Nf], sqrt(diag(UH[:Nf,:Nf])),fmt=".-",color=colors[2],alpha=0.2)
xlabel("frequency / MHz",fontsize=20); ylabel("real part / au",fontsize=20)
subplot(122)
errorbar(f*1e-6, H[Nf:],sqrt(diag(UH[Nf:,Nf:])),fmt=".-",color=colors[2],alpha=0.2)
xlabel("frequency / MHz",fontsize=20); ylabel("imaginary part / au",fontsize=20);
The deconvolution problem can be decomposed into a division by the system's frequency response followed by a multiplication by a low-pass filter frequency response.
X(f) = \frac{Y(f)}{H(f)}H_L(f)
which in real and imaginary part becomes
X = \frac{(\Re_Y\Re_H + \Im_Y\Im_H)+\mathrm{j}(-\Re_Y\Im_H+\Im_Y\Re_H)}{\Re_H^2+\Im_H^2}(\Re_{H_L} + \mathrm{j}\Im_{H_L})
Sensitivities for division part
Uncertainty blocks for multiplication part
# low-pass filter for deconvolution
def lowpass(f,fcut=80e6):
return 1/(1+1j*f/fcut)**2
HLc = lowpass(f)
HL = np.r_[np.real(HLc), np.imag(HLc)]
XH,UXH = DFT_deconv(H,Y,UH,UY)
XH, UXH = DFT_multiply(XH, UXH, HL)
figure(figsize=(18,6))
subplot(121)
errorbar(f*1e-6, XH[:Nf], sqrt(diag(UXH[:Nf,:Nf])),fmt=".-",color=colors[2],alpha=0.2)
xlabel("frequency / MHz",fontsize=20); ylabel("real part / au",fontsize=20)
subplot(122)
errorbar(f*1e-6, XH[Nf:],sqrt(diag(UXH[Nf:,Nf:])),fmt=".-",color=colors[2],alpha=0.2)
xlabel("frequency / MHz",fontsize=20); ylabel("imaginary part / au",fontsize=20);
The inverse DFT equation is given by
X_n = \frac{1}{N} \sum_{k=0}^{N-1}(\Re_k\cos(k\beta_n)-\Im_k\sin(k\beta_n))
The sensitivities for the GUM propagation of uncertainties are then
GUM uncertainty propagation for the inverse DFT
xh,Uxh = GUM_iDFT(XH,UXH,Nx=N)
ux = np.sqrt(np.diag(Uxh))
figure(figsize=(16,8))
plot(time*1e6,xh,label="estimated pressure signal",linewidth=2,color=colors[0])
plot(time * 1e6, ref_data, "--", label= "reference data", linewidth=2, color=colors[1])
fill_between(time * 1e6, xh + 2 * ux, xh - 2 * ux, alpha=0.2, color=colors[0])
xlabel("time / µs", fontsize=22)
ylabel("signal amplitude / MPa", fontsize=22)
tick_params(which= "major", labelsize=18)
legend(loc= "upper left", fontsize=18, fancybox=True)
xlim(0, 2);
figure(figsize=(16,8))
plot(time * 1e6, ux, label= "uncertainty", linewidth=2, color=colors[0])
xlabel("time / µs", fontsize=22)
ylabel("uncertainty / MPa", fontsize=22)
tick_params(which= "major", labelsize=18)
xlim(0,2);
Y,UY = GUM_DFT(y,Uy,N=Nf) H, UH = AmpPhase2DFT(A, P, UAP) XH,UXH = DFT_deconv(H,Y,UH,UY) XH, UXH = DFT_multiply(XH, UXH, HL)